package scu.maqiang.electromagnetics.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarCFunc;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Complex;
import scu.maqiang.numeric.IterSCSolver;
import scu.maqiang.numeric.SCMatrix;
import scu.maqiang.numeric.SRMatrix;

public class Scattering1D_Fine {
    public static void main(String[] args) {
        double c0 = 3.0e8; //自由空间光速
        double nu0 = 120 * Math.PI; //自由空间固有阻抗
        double e0 = (1e-9)/(36*Math.PI); //自由空间电容率
        double mu0 = 4 * Math.PI * 1e-7; //自由空间在磁导率

        double freq = 300 * 1e6; //频率
        double lambda = c0 / freq; //自由空间波长
        double k0 = 2 * Math.PI / lambda; //自由空间波数
        double omg = 2 * Math.PI * freq; //径向频率

        double er = 2.7; //板电容率
        double sigma = 5e-3; //板电导率
        double mur = 1;  //板相对磁导率(非磁体介质)

        double L = 2.0;    //区域长度
        double d = 1.0;   //板厚度
        double xa = 0.5 * (L - d);  //板左端坐标
        double xb = 0.5 * (L + d);  //板右端坐标
        Mesh1L2 mesh = new Mesh1L2();
        mesh.line(0, L, 160);
        ScalarFunc slabRegion = (xy, label, param) -> (xy[0] > xa && xy[0] < xb)? 1.0: -1.0;
        int slabRegionNumber = 10;
        mesh.setDomainLabel(slabRegion, null, slabRegionNumber);

        Complex erc = new Complex(-k0 * k0 * er, k0 * k0 * sigma /(omg * e0));
        System.out.println(erc);
        Complex er0 = new Complex(-k0 * k0, 0.0);
        FES1L21 fs = new FES1L21(mesh);
        ScalarFunc muCoef = (xy, label, param) -> label == slabRegionNumber? 1.0 / mur: 1.0;
        SRMatrix A = new SRMatrix(fs.GetNdof());
        fs.assembleStiff(new ScalarFunc[]{muCoef}, null, BVPType.COMMON, A);

        ScalarCFunc eCoef = (xy, label, param) -> label == slabRegionNumber? erc: er0;
        SCMatrix M = new SCMatrix(fs.GetNdof());
        fs.assembleCMass(new ScalarCFunc[]{eCoef}, null, BVPType.CONSIST_MASS, M);

        SCMatrix K = M.add(1.0, A);
//        K.sort();
//        System.out.println(K.display());

        fs.applyABC(K, k0 / 1.0, 1, 2);
        //System.out.println(K.display());

        ScalarCFunc sourceCoef = (xy, label, param) -> {
            double realpart = k0 * k0 * ((er - 1) * Math.cos(k0 * xy[0]) - sigma / (omg * e0) * Math.sin(k0 * xy[0]));
            double imagpart = -k0 * k0 * ((er - 1) * Math.sin(k0 * xy[0]) + sigma / (omg * e0) * Math.cos(k0 * xy[0]));
            return label == slabRegionNumber? new Complex(realpart, imagpart):new Complex(0, 0);
        };

        //        double[][]  xy = new double[2][2];
//        for(int i = 0; i < 80; i++) {
//            mesh.getCoordPerElement(mesh.getElement(i), xy);
//            double x = 0.5 * (xy[0][0] + xy[1][0]);
//            System.out.print("-------------" + i + "\t   ");
//            System.out.print(Math.cos(k0 * x) + "\t\t" + Math.sin(k0 * x) + "\t  ");
//            System.out.print(k0 * k0 * ((er - 1) * Math.cos(k0 * x) - sigma / (omg * e0) * Math.sin(k0 * x)) + "\t   ");
//            System.out.println(-k0 * k0 * ((er - 1) * Math.sin(k0 * x) +sigma / (omg * e0) * Math.cos(k0 * x)));
//        }
        Complex[] RHS = Complex.zeros(fs.GetNdof());
        fs.assembleCSource(new ScalarCFunc[]{sourceCoef}, null, BVPType.COMMON, RHS);
        System.out.println(Complex.toString(RHS));

        IterSCSolver solver = new IterSCSolver(K);
        Complex[] u = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR(RHS, u, 1.5, 1);


        double[][] result = new double[4][];
        result[0] = Complex.Real(u);
        result[1] = Complex.Image(u);
        result[2] = Complex.abs(u);

        ScalarCFunc incidentField = (xy, label, param) -> new Complex(Math.cos(k0 * xy[0]), -Math.sin(k0 * xy[0]));
        Complex[] uincn = mesh.functionInMesh(incidentField, null);
        Complex[] utot = Complex.add(uincn, 1.0, u);
        result[3] = Complex.abs(utot);
        mesh.toTecplot("Scattering1D_Fine.dat", result);
    }

}
